

/*
 * Class:        CIRProcess
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */

package umontreal.iro.lecuyer.stochprocess;
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;



/**
 * This class represents a <SPAN  CLASS="textit">CIR</SPAN> (Cox, Ingersoll, Ross) process
 *  
 * <SPAN CLASS="MATH">{<I>X</I>(<I>t</I>) : <I>t</I>&nbsp;&gt;=&nbsp;0}</SPAN>, sampled at times 
 * <SPAN CLASS="MATH">0 = <I>t</I><SUB>0</SUB> &lt; <I>t</I><SUB>1</SUB> &lt; <SUP> ... </SUP> &lt; <I>t</I><SUB>d</SUB></SPAN>.
 * This process obeys the stochastic differential equation
 * 
 * <P></P>
 * <DIV ALIGN="CENTER" CLASS="mathdisplay"><A NAME="eq:cir"></A>
 * <I>dX</I>(<I>t</I>) = <I>&#945;</I>(<I>b</I> - <I>X</I>(<I>t</I>))<I>dt</I> + <I>&#963;</I>(X(t))<SUP>1/2</SUP>&nbsp;<I>dB</I>(<I>t</I>)
 * </DIV><P></P>
 * with initial condition <SPAN CLASS="MATH"><I>X</I>(0) = <I>x</I><SUB>0</SUB></SPAN>,
 * where <SPAN CLASS="MATH"><I>&#945;</I></SPAN>, <SPAN CLASS="MATH"><I>b</I></SPAN> and <SPAN CLASS="MATH"><I>&#963;</I></SPAN> are positive constants,
 * and 
 * <SPAN CLASS="MATH">{<I>B</I>(<I>t</I>),&nbsp;<I>t</I>&nbsp;&gt;=&nbsp;0}</SPAN> is a standard Brownian motion
 * (with drift 0 and volatility 1).
 * This process is <SPAN  CLASS="textit">mean-reverting</SPAN> in the sense that it always tends to
 * drift toward its general mean <SPAN CLASS="MATH"><I>b</I></SPAN>.
 * The process is generated using the sequential
 * technique
 * 
 * <P></P>
 * <DIV ALIGN="CENTER" CLASS="mathdisplay"><A NAME="eq:cir-seq"></A>
 * <I>X</I>(<I>t</I><SUB>j</SUB>) = {<I>&#963;</I><SUP>2</SUP>(1-<I>e</I><SUP>-<I>&#945;</I>(t<SUB>j</SUB>-t<SUB>j-1</SUB>)</SUP>)/4<I>&#945;</I>}<I>&#967;</I><SUP>&#8242;2</SUP><SUB><I>&#957;</I></SUB>({4<I>&#945;e</I><SUP>-<I>&#945;</I>(t<SUB>j</SUB>-t<SUB>j-1</SUB>)</SUP><I>X</I>(<I>t</I><SUB>j-1</SUB>)}/{<I>&#963;</I><SUP>2</SUP>(1-<I>e</I><SUP>-<I>&#945;</I>(t<SUB>j</SUB>-t<SUB>j-1</SUB>)</SUP>)}),
 * </DIV><P></P>
 * where 
 * <SPAN CLASS="MATH"><I>&#957;</I> = 4<I>b&#945;</I>/<I>&#963;</I><SUP>2</SUP></SPAN>, and 
 * <SPAN CLASS="MATH"><I>&#967;</I><SUP>&#8242;&nbsp;2</SUP><SUB><I>&#957;</I></SUB>(<I>&#955;</I>)</SPAN> is a noncentral
 * chi-square random variable with <SPAN CLASS="MATH"><I>&#957;</I></SPAN> degrees of freedom and noncentrality
 * parameter <SPAN CLASS="MATH"><I>&#955;</I></SPAN>.
 * 
 */
public class CIRProcess extends StochasticProcess  {
    private RandomStream stream;
    protected ChiSquareNoncentralGen gen;
    protected double alpha,
                     beta,
                     sigma,
                     nu;     // number of degrees of freedom
    // Precomputed values
    protected double[] parc,
                       parlam;



   /**
    * Constructs a new <TT>CIRProcess</TT> with parameters
    * 
    * <SPAN CLASS="MATH"><I>&#945;</I> = <texttt>alpha</texttt></SPAN>, <SPAN CLASS="MATH"><I>b</I></SPAN>, 
    * <SPAN CLASS="MATH"><I>&#963;</I> = <texttt>sigma</texttt></SPAN> and initial value
    * 
    * <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>0</SUB>) = <texttt>x0</texttt></SPAN>. The noncentral chi-square variates
    *   
    * <SPAN CLASS="MATH"><I>&#967;</I><SUP>&#8242;2</SUP><SUB><I>&#957;</I></SUB>(<I>&#955;</I>)</SPAN> will be
    * generated by inversion using the stream <TT>stream</TT>.
    * 
    */
   public CIRProcess (double x0, double alpha, double b, double sigma,
                      RandomStream stream)  {
      this (x0, alpha, b, sigma, new ChiSquareNoncentralGen (stream, null));
    }


   /**
    * The noncentral chi-square variate generator  <TT>gen</TT>
    *  is specified directly instead of specifying the stream.
    *  <TT>gen</TT> can use a method other than inversion.
    * 
    */
   public CIRProcess (double x0, double alpha, double b, double sigma,
                      ChiSquareNoncentralGen gen)  {
      this.alpha = alpha;
      this.beta  = b;
      this.sigma = sigma;
      this.x0    = x0;
      nu = 4.0*b*alpha/(sigma*sigma);
      this.gen   = gen;
      stream = gen.getStream();
    }


   public double nextObservation() {
      double xOld = path[observationIndex];
      double lambda = xOld * parlam[observationIndex];
      double x;
      if (gen.getClass() == ChiSquareNoncentralPoisGen.class)
         x = parc[observationIndex] *
             ChiSquareNoncentralPoisGen.nextDouble(stream, nu, lambda);
      else if (gen.getClass() == ChiSquareNoncentralGamGen.class)
         x = parc[observationIndex] *
             ChiSquareNoncentralGamGen.nextDouble(stream, nu, lambda);
      else
         x = parc[observationIndex] *
             ChiSquareNoncentralGen.nextDouble(stream, nu, lambda);
      observationIndex++;
      path[observationIndex] = x;
      return x;
   }


   /**
    * Generates and returns the next observation at time 
    * <SPAN CLASS="MATH"><I>t</I><SUB>j+1</SUB> = <texttt>nextTime</texttt></SPAN>, using the previous observation time <SPAN CLASS="MATH"><I>t</I><SUB>j</SUB></SPAN> defined earlier
    * (either by this method or by <TT>setObservationTimes</TT>),
    * as well as the value of the previous observation <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>j</SUB>)</SPAN>.
    * <SPAN  CLASS="textit">Warning</SPAN>: This method will reset the observations time <SPAN CLASS="MATH"><I>t</I><SUB>j+1</SUB></SPAN>
    * for this process to <TT>nextTime</TT>. The user must make sure that
    * the <SPAN CLASS="MATH"><I>t</I><SUB>j+1</SUB></SPAN> supplied is 
    * <SPAN CLASS="MATH">&nbsp;&gt;=&nbsp;<I>t</I><SUB>j</SUB></SPAN>.
    * 
    */
   public double nextObservation (double nextTime)  {
      double previousTime = t[observationIndex];
      double xOld = path[observationIndex];
      observationIndex++;
      t[observationIndex] = nextTime;
      double dt = nextTime - previousTime;
      double x = nextObservation (xOld, dt);
      path[observationIndex] = x;
      return x;
    }


   /**
    * Generates an observation of the process in <TT>dt</TT> time units,
    * assuming that the process has value <SPAN CLASS="MATH"><I>x</I></SPAN> at the current time.
    * Uses the process parameters specified in the constructor.
    * Note that this method does not affect the sample path of the process
    * stored internally (if any).
    * 
    * 
    */
   public double nextObservation (double x, double dt)  {
      double c = -Math.expm1(-alpha * dt) * sigma * sigma / (4.0*alpha);
      double lambda = x * Math.exp(-alpha * dt) / c;
      if (gen.getClass() == ChiSquareNoncentralPoisGen.class)
         x = c * ChiSquareNoncentralPoisGen.nextDouble(stream, nu, lambda);
      else if (gen.getClass() == ChiSquareNoncentralGamGen.class)
         x = c * ChiSquareNoncentralGamGen.nextDouble(stream, nu, lambda);
      else
         x = c * ChiSquareNoncentralGen.nextDouble(stream, nu, lambda);
      return x;
    }


   public double[] generatePath() {
      double xOld = x0;
      double x, lambda;
      int j;

      if (gen.getClass() == ChiSquareNoncentralPoisGen.class) {
         for (j = 0; j < d; j++) {
            lambda = xOld * parlam[j];
            x = parc[j] * ChiSquareNoncentralPoisGen.nextDouble(stream, nu, lambda);
            path[j + 1] = x;
            xOld = x;
         }

      } else if (gen.getClass() == ChiSquareNoncentralGamGen.class) {
         for (j = 0; j < d; j++) {
            lambda = xOld * parlam[j];
            x = parc[j] * ChiSquareNoncentralGamGen.nextDouble(stream, nu, lambda);
            path[j + 1] = x;
            xOld = x;
         }

      } else {
         for (j = 0; j < d; j++) {
            lambda = xOld * parlam[j];
            x = parc[j] * ChiSquareNoncentralGen.nextDouble(stream, nu, lambda);
            path[j + 1] = x;
            xOld = x;
         }
      }

      observationIndex = d;
      return path;
   }

   public double[] generatePath (RandomStream stream) {
      gen.setStream (stream);
      return generatePath();
   }


   /**
    * Resets the parameters 
    * <SPAN CLASS="MATH"><I>X</I>(<I>t</I><SUB>0</SUB>) = <texttt>x0</texttt></SPAN>, 
    * <SPAN CLASS="MATH"><I>&#945;</I> = <texttt>alpha</texttt></SPAN>,
    *  
    * <SPAN CLASS="MATH"><I>b</I> = <texttt>b</texttt></SPAN> and 
    * <SPAN CLASS="MATH"><I>&#963;</I> = <texttt>sigma</texttt></SPAN> of the process.
    * <SPAN  CLASS="textit">Warning</SPAN>: This method will recompute some quantities stored internally,
    * which may be slow if called too frequently.
    * 
    */
   public void setParams (double x0, double alpha, double b, double sigma)  {
      this.alpha = alpha;
      this.beta  = b;
      this.sigma = sigma;
      this.x0    = x0;
      nu = 4.0*b*alpha/(sigma*sigma);
      if (observationTimesSet) init(); // Otherwise not needed.
    }


   /**
    * Resets the random stream of the noncentral chi-square generator to <TT>stream</TT>.
    * 
    */
   public void setStream (RandomStream stream)  { gen.setStream (stream); }


   /**
    * Returns the random stream of the noncentral chi-square generator.
    * 
    */
   public RandomStream getStream()  { return gen.getStream (); }


   /**
    * Returns the value of <SPAN CLASS="MATH"><I>&#945;</I></SPAN>.
    * 
    */
   public double getAlpha()  { return alpha; }


   /**
    * Returns the value of <SPAN CLASS="MATH"><I>b</I></SPAN>.
    * 
    */
   public double getB()  { return beta; }


   /**
    * Returns the value of <SPAN CLASS="MATH"><I>&#963;</I></SPAN>.
    * 
    */
   public double getSigma()  { return sigma; }


   /**
    * Returns the noncentral chi-square random variate generator used.
    * The <TT>RandomStream</TT> used for that generator can be changed via
    * <TT>getGen().setStream(stream)</TT>, for example.
    * 
    */
   public ChiSquareNoncentralGen getGen()  { return gen; }
 

   protected void initArrays(int d) {
      double dt, c;
      for (int j = 0; j < d; j++) {
         dt = t[j+1] - t[j];
         c = -Math.expm1(-alpha * dt)*sigma*sigma/(4.0*alpha);
         parc[j] = c;
         parlam[j] = Math.exp(-alpha * dt) / c;
      }
   }

   // This is called by setObservationTimes to precompute constants
   // in order to speed up the path generation.
   protected void init() {
       super.init();
       parc = new double[d];
       parlam = new double[d];
       initArrays(d);
    }

} 
